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Based on the analytic model of Feshbach resonances in harmonic traps described in Phys. Rev. A 
83, 030701 (2011) a Bose-Hubbard model is introduced that provides an accurate description of two 
atoms in an optical lattice at a Feshbach resonance with only a small number of Bloch bands. The 
approach circumvents the problem that the eigenenergies in the presence of a delta-like coupling do 
not converge to the correct energies, if an uncorrelated basis is used. The predictions of the Bose- 
Hubbard model are compared to non-perturbative calculations for both the stationary states and 
the time-dependent wavefunction during an acceleration of the lattice potential. For this purpose, 
a square- well interaction potential is introduced, which allows for a realistic description of Feshbach 
resonances within non-perturbative single-channel calculations. 



I. INTRODUCTION 

Since the creation of the first Bosc-Einstein conden- 
sates [H [2] , ultracold atoms have proven to be a versatile 
tool for many applications like precision measurement, 
quantum simulation, and quantum information process- 
ing. Two of the main techniques that made these achieve- 
ments possible are the creation of various trapping po- 
tentials, like optical lattices (OLs) or wave guides and the 
precise control of the interatomic interaction by means of 
Feshbach resonances (FRs) 0] . 

An important tool for describing ultracold atoms in 
OLs is the Bose-Hubbard (BH) model. The model uses 
in its basic form a basis of single-particle Wannier states 
from the first Bloch band to formulate the many-body 
Hamiltonian. While for weak interactions the model is 
very accurate, it usually breaks down for larger scattering 
lengths. A way to extend its applicability at a broad FR 
is to introduce effective BH parameters especially for the 
onsite interaction strength U. These parameters can be 
obtained by using a corrected harmonic approximation 
of the lattice sites [5] or by full numerical calculations 

mm- 

The usual BH model allows via the onsite-interaction 
strength U either for repulsively interacting atoms (U > 
0) or attractively interacting atoms (U < 0). At a narrow 
FR, however, a relatively narrow avoided crossing with 
the resonant bound state leads to the appearance of both 
repulsively and attractively interacting states [HI IS] ■ In 
this situation the resonant bound state must be explicitly 
included into the BH model. Several different kinds of 
these extended models have been introduced and debated 
[5HT2] and applied to map out the phase diagramm [TU1 
fl3l IT4] or to investigate lattice solitons [15] . 

The above investigations consider the extended Hub- 
bard model within a single-band approximation that is 
only applicable in the rare situation that the coupling 
energy to the resonant bound state is small compared to 
the band gap. In order to generalize the applicability 
one can introduce the notion of dressed molecules with 
effective bound-state energies and coupling strengths ob- 
tained from more elaborate calculations [TB] . 

A convenient approach to generalize Hubbard models 



to describe broader FRs or systems with a large scat- 
tering length is to simply include more Bloch bands. For 
example, L.-M. Duan has derived an effective single-band 
Hubbard model for the case of interacting fermions at a 
broad FR starting from a multi-band Hubbard model 
in the Wannier basis and a zero-range coupling between 
atoms and molecules [T7J. However, as will be discussed 
in this work, severe numerical problems arise for the de- 
scription of a system with a zero-range coupling, e.g., by 
expanding the solution in products of single-particle ba- 
sis functions. Especially for large scattering lengths all 
of these basis functions behave completely differently for 
r — > compared with the correct solution. This poses a 
problem especially for positive scattering lengths where 
the open channel supports a bound state. In fact, the 
obtained energies are lower than the correct ones so that 
an increase of the basis leads to an even larger disagree- 
ment. A similar problem also appears when replacing 
the interaction potential by the delta-like Fermi-Huang 
pseudo-potential |18j . Also within analytical treatments 
of FRs in harmonic traps that use non-interacting basis 
states the eigenenergies do not converge jTDl [TH] ■ In this 
case, after an infinite summation, the diverging terms 
can be absorbed by introducing a renormalized bound- 
state energy. In many numerical approaches the problem 
is circumvented by replacing the delta-like potential by 
a regularized short range potential [H [7J [5D] . In order 
to resolve the potential usually a large basis is neces- 
sary. For example, for an interaction with the range d/N 
where d is the lattice spacing more than N Bloch bands 
have to be included to converge the energies [3 [7J . Since 
for two atoms in a one dimensional lattice the number 
of basis functions scales quadratically with the number 
of Bloch bands and the number of sites the solution can 
quickly become numerically very demanding. Based on 
this corrected numerical approach, M.L. Wall and L.D. 
Carr were able to calculated the effective parameters of 
a Fermi Hubbard model that takes the coupling to a 
bosonic molecule explicitly into account (16) . 

In this work we introduce an extended BH model that 
avoids the numerical problems in the presence of a delta- 
like coupling without the need of regularization and in- 
clusion of many Bloch bands. The model is derived from 



2 



first principles on the basis of the analytic microscopic 
theory of FRs in a harmonic trap [8]. This allows for 
defining dressed bound-state energies and couplings that 
correct for the problems due to the deficiency of the basis 
states. 

Given the number of different proposals to describe 
FRs within a BH model one has to compare the predic- 
tions of the introduced BH model with non-perturbative 
calculations. In the standard description of FRs this re- 
quires to solve a two-channel problem of two interacting 
atoms in an optical lattice coupled at short distance to 
a molecular bound state. This problem is numerically 
very demanding. However, we show that one can largely 
simplify the problem by introducing a square-well inter- 
action potential that realistically mimics the behavior at 
a FR. Using this single-channel interaction potential we 
apply an approach introduced in [211 122) in order to ob- 
tain the correct energies and wave functions of two atoms 
in a small OL at a FR. The correct stationary and dy- 
namic behavior of two atoms in a double-well potential is 
compared with the results of the introduced BH model. 
It is shown that with only a small number of Bloch bands 
included the BH model is able to accurately describe FRs 
of small and medium width with coupling energies up to 
the depth of the OL. 

The work is organized as follows. First the analytic 
model of a FR in a harmonic introduced in [5] is briefly 
recapitulated. This sets the basis for the derivation of a 
general BH model of interacting atoms at a FR in Sec. |III| 
The model is compared to the exact analytical solution in 
a harmonic trap in Sec. |IV| revealing that the BH model 
does not converge toward the correct eigenenergies. To 
circumvent this problem dressed molecular states and a 
dressed coupling strength of the BH model are introduced 
in Sec. [V] on the basis of the analytically known eigenen- 
ergies in the harmonic trap. In Sec. |VI| the square-well 
interaction potential is discussed, which allows for finding 
within a non-perturbative approach both the stationary 
and the time-dependent wavefunctions of two atoms in a 
small OL at a FR perturbed by a time-dependent accel- 
eration of the lattice. Finally, in Sec. |VII| the dressed and 
undressed BH model is compared to the non-perturbative 
calculations. We conclude in Sec. I Villi 



II. FESHBACH RESONANCE IN A HARMONIC 
TRAP 

Neutral atoms usually only interact at small distances 
r on the order of r; nt ~ 100 a.u. which is much smaller 
than typical length scales of the trapping potentials on 
the order of some r trap ~ 10 000 a.u. The collision energy 
in the ultracold regime is so small that partial waves with 
angular momentum I > arc reflected by the centrifugal 
barrier. Therefore s-wave (I — 0) scattering is largely 
dominant. For r mt <C r <C r tra p the interaction leads to 
a phase shift ip of the scattering wave function oc sin(fcr + 
ip) which is associated with the s-wave scattering length 



a(k) = — tan(<^)/fc. 

If two atoms collide, the spin states of the scattering 
atoms are coupled at small distances r < r; nt to other 
spin states in closed channels whose relative energy can 
be influenced by applying an external magnetic field B. 
The subspace of closed-channel spin states can support 
many bound states. For certain magnetic field strengths 
B the energy of such a bound state E^B) can be brought 
into resonance with the collision energy E of the atoms, 
leading to a FR (see Fig. [I]). 
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FIG. 1. (Color online). Mechanism of broad and narrow 
FRs in the exemplary case of two atoms in isotropic har- 
monic confinement with frequency u. Top: Sketch of the 
relative-motion unbound trap states and the resonant molec- 
ular bound state (RBS) whose energy can be manipulated by 
an external magnetic field. Bottom: Relative-motion energy 



spectrum [solutions of Eq. 1 12 1] using the energy-dependent 
scattering length a(E, E TeB ) (thick red lines) and the energy- 
independent scattering length a(0, -E res ) (black dashed line) 
as a function of the resonance energy E ICS — Et + 8E (black 
solid line). At a narrow FR (left: ab g = 0.04ah o , A_E = lhco) 
the RBS couples only to the trap state that is in resonance, 
which leads to narrow avoided crossings. At a broad FR 
(right: a bg = 0.04a ho , A_E = 40huj) the RBS couples to many 
trap states and the energy spectrum changes globally with 
E ICS . In contrast to the narrow FR the eigenenergies for an 
energy-dependent and energy-independent scattering length 
agree reasonably for broad resonances. 
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In [5] an analytic model for a FR in isotropic and 
anisotropic harmonic traps was developed. Its starting 
point is the relative-motion (REL) Hamiltonian for radial 
momentum I = of two atoms in a spherical harmonic 
confinement with frequency u). The Hamiltonian for the 
radial wave function = V Anr |$(r)), where \^(r)) 

is the REL wave function, is given as 



H 



K 1 d 2 



2^ + 2^ V+VzH+V - W - (1) 



Here, p is the reduced mass, Vzh is the Zeeman and 
hyperfine energy of the atoms, and Vi nt (r) the electron- 
spin dependent interaction potential. 

One assumes that the REL energy of the atoms is small 
enough so that only one spin configuration (the open 
channel) supports unbound states. All other spin con- 
figurations are closed, i.e. their wave function vanishes 
for large interatomic distances. 

The analytic model is based on the two-channel de- 
scription of an FR where one channel represents the 
unbound atoms and the other channel the subspace of 
closed channels. By introducing projectors P and Q onto 
the subspace of open and closed channels, respectively, 
one arrives at the coupled equations 



(Hp - E) \<h P ) 
(fig - E) |$q> 



fW|$g) 
W f |$ P ) 



(2) 
(3) 



with Hp = PHP, H Q = QHQ, W = PHQ, |$ P ) = P |$), 
I^q) = Ql^)i and E the energy above the thresh- 
old of the open-channel interaction potential. Since the 
eigenenergies in the closed channel subspace are usually 
largely separated on the energy scale hui of the trap, one 
assumes that close to the FR |$q) is simply a multiple 
A of a single bound eigenstate \<pb) with eigenenergy E b . 
We call this closed-channel state "resonant bound state" 
(RBS). To first order, the energy E b may be expanded 
linearly in the magnetic field B, i.e. Eb(B) = a(B — Bo), 
where a is the relative magnetic moment that is known 
for many FRs [4]. 

Introducing the normalized solution \4>e) of the open 
channel with \<&p) — C \4>e) and a background eigenstate 
1 0bg) of Hp with eigenenergy E\> g which is occupied for 
infinite detuning \E — Eb\ — > oo one obtains the eigenen- 
ergy equation [23] 



(E - E b )(E - E bg ) = 



W| 



iwt 



9 bg| <PE) 



(4) 



In order to find simplified expressions for (</>bg| W \ 4>b), 
(<pb\ \4>e)i and (</>b g | 4*e) one assumes that the inter- 
action is only relevant in some small range r < rj nt much 
smaller than the extension of the trap r tra p. The exten- 
sion of the harmonic trap is specified by the harmonic 
trap length a^o = yV (pui). Denoting the long-range 
behavior of ^js(r) by <pE{r) one finds 



4>s(r) = lim <j> E {r) = A v D v (p), 



(5) 



where D v (p) is the parabolic cylinder function, p = 
v2r/aho! v — E / (Fkj) — 1/2, and A v is a normalization 
constant. 

For r <C a\ w the linear approximation of D v (p) 
yields [H] 



MO) 1 



with 



f(E) 



a ho f(E) 

1 V4 2hu) 

or (2 _ 



+ 0(r 2 ) 



(6) 



(7) 



where F(x) is the Gamma function. In the range r int <C 
r <C aho the radial wavefunction with scattering length a 
has the form 4>E{r) cx 1—r/a. Hence, one can directly de- 
termine the scattering length of the radial wave function 
with energy E from Eq. Q. This yields 



;(0) 



= a ho f(E) 



(8) 



which is equivalent to the result in [25] , 

In the spirit of a Taylor expansion we parametrize 
(<Pe\ W \4>b) by a linear combination 



(0 6 |Wt|<M=7MO)+/^ / E (O) 
= 7<M0) + 



^(0) 



7 MO)/ ( 9 ) 



7<M0) 1- 



with a* — /3/j. Be ipb(r) = 4>b { r )/{ v47rr) the wave 
function describing the RBS then the expansion ^ can 
be interpreted as approximating the coupling to the 
bound state by W(r)ipb(r) ~ \^4n^(r — a*)5(f). For 
the long-range behavior of the wavefunction ifjE(r), i.e. 
linV-^oo ipE^r) = rp E {r) = 4>E(r)/{V^r) one finds 

7^(0)+/3<^(0) 

dr4> E (r) [*/8(r)+j38'(r)] 

(10) 

dry 47T'0p(r) [77™ — /3] S(r) 



r dr dfl ipE( r )y 4tt7 [r — a*] S(r) 



Here one uses rS'(r) = —S(r) and <5(r) = 4irr 2 S(r). Al- 
though only two parameters are used, the parametriza- 
tion of the coupling is already quite general since higher 
order couplings like those proportional to r 2 5(r) auto- 
matically vanish. Within the approximation of a con- 
stant RBS 7 and o* must be constant. In reality, how- 
ever, (4>b I \4>e) depends on the nodal structure of the 
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RBS and the open channel that are both not constant 
for a varying magnetic field. A comparison with com- 
plete coupled-channel calculations shows that it suffices 
to introduce a background coupling strength 7b g for the 

parametrization of (4>bg\ W \ 4>b) to account for slight vari- 
ations of the nodal structure [5] . Since the difference be- 
tween 7 and 7b g is only relevant for the RBS admixture 
but not for the eigenenergies of the system, we can ignore 
it for our purposes. Following the reasoning given in [8], 
the short-range approximation ^ then gives 



E-E h = 



2 7 2 {fw £) (/(^) - £) 



f(E) - f(E hg ) 



(11) 



The solutions of this equation determine the eigenener- 
gies. One can rewrite this equation in the form of a 
matching condition: The scattering length a(E, E^) due 
to the short-range coupling to the RBS must be equal 
to the product a\ lo f(E) that is equal to the scattering 
length of the long-range wavefunction </>s(r). This yields 



a ho f(E) = a(E,E 

res) — G&bg iE) 1 



AE 



E IC 



(12) 



The right-hand side of the Eq. ( 12 ) describes the energy 



dependence of the scattering length with background 
scattering length ab g and resonance width 



h 2 V a bg 



(13) 



The resonance energy E Ies = E^ + SE is shifted from the 
bound state energy Ef, by the resonance detuning 



5E = 



a\,„AE 



a bg 



a* 



(14) 



In the limit E — > the ratio of the resonance de- 
tuning and the resonance width is given as SE/AE = 
ao/(io — where ao = lim^^oabg is the zero-energy 
background scattering length. Comparing this with the 
same ratio derived on the basis of a multi-channel quan- 
tum defect theory for E — > 26 allows for removing one 
free parameter a*. One finds 



1 



ao 



(15) 



where the mean scattering length a is determined by the 
Cq coefficient of the van der Waals interaction [57] . Using 
Eqs. (13) and (|15[), the remaining parameter 7 can be 



directly related to the resonance width AE. 

The function a- ilo f{E) which describes the scattering 
length of the wave function ^(r) is also known for 
anisotropic traps with u> y — u z = 7]lj x . In this case the 
scattering length is given as a = —y/Trd/F(u,r]) (with 
d, u = u(E), and T defined in [28]) such that the eigenen- 
ergy relation 



fixd 



= a hs (E) 1- 



AE 



E IC 



E 



(16) 



holds. 

One generally distinguishes between narrow and broad 
FRs [1] . In the case of a broad FR the coupling strength 
to the bound state is relatively large such that it is 
admixed to unbound states in a large energy domain. 
If, as usual, the background scattering length a bg is 
small compared to the trap length ah , the ratio of the 
RBS admixture \A\ 2 to the open-channel admixture |C| 2 
for states above the first trap state is on the order of 
CLhoHw/ (ahgAE) such that the RBS admixture can be 
neglected if a^ g AE S> a\ m huj 8 . Furthermore, also the 
energy dependence of the scattering length becomes neg- 
ligible if a,b g AE 3> a\i huj (see Fig. [IJ. Therefore, all 
details of the atomic interaction apart from the value of 
scattering length for E — > are irrelevant. This situation 
is called universal. On the other hand, for narrow FRs 
the bound state couples only to a narrow energy range 
of scattering states or respectively to that unbound trap 
state which is in resonance. As shown in Fig. [T]in a har- 
monic trap this leads to narrow avoided crossings in the 
energy spectrum with an energy splitting on the order 
of ^a hg AE I (a ho huj)hu) [5]. At the resonance the bound 
state is strongly admixed and the energy dependence of 
the scattering length cannot be neglected. 



III. FESHBACH RESONANCE IN AN OPTICAL 
LATTICE 

In order to avoid unnecessary complexity, in the follow- 
ing an OL is considered, in which two directions of move- 
ment are effectively frozen out by using strong harmonic 
confinement. Nevertheless, the following discussions can 
be easily extended to 2D and 3D lattices. 

A particle of mass m in such an OL of depth Vl and pe- 
riodicity d = ir/ko in the spacial direction x and transver- 
sal harmonic confinement with frequency uj t in y and z 
direction is described by the Hamiltonian 



P 

2m 



+ V L sin 2 (fc i) + xmw t 2 (y 2 + z 2 ) . (17) 



Eigensolutions of this Hamiltonian with quasi momentum 
k are given by 

®k,n,M v ,Mz ( x > V' z ) = e ikx <t>n,k{x) H My (y) H Mz (z) (18) 

where <p ni k are analytically known Bloch solutions with 
band index n — 1,2,3,... and quasi momentum k of 
the periodic lattice. Hm is the M-th solution of the 
one-dimensional harmonic oscillator with transversal fre- 
quency Wf 

In order to describe more than one particle in an OL, 
interactions have to be taken into account. Since neutral 
atoms interact only on short distances it is convenient to 



transform the basis ( 18 ) into localized functions. This is 



done by the usual transformation to Wannier functions 



Wi <n< M v , M Ax,V,z) = Wi,n{x)H M Jy)H Mx {z). (19) 
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Here, Wi jn denotes the Wannier function localized at lat- 
tice site i and band n. 

Due to the anharmonicity of the OL the relative- 
motion (REL) coordinates f — (x, y, z) T = r\ — r 2 and 
center-of-mass (COM) coordinates R — (X,Y,Z) T = 
(ri +r2)/2 are coupled. Therefore, the Eqs. ([2| and ([3| 
for REL motion have to be extended to include also the 
COM energies of the two atoms and the resonant molec- 
ular state. To this end &p(ri, f 2 ) shall describe the wave 
function of the two atoms in the open channel with ki- 
netic and potential energies / H m (fi) +W m ('r f 2) interacting 
via a short-range potential V(r). The open channel is 
coupled by some real- valued short-range coupling W(r) 
to the closed-channel wave function ^q(R, r). One as- 
sumes that the RBS in REL motion has an extension 
small enough not to probe the external trapping poten- 
tial. Therefore, the closed-channel wave function can be 
written as a product state ^q(R, r) = V'bfT )^COm(-R) of 
the RBS ipb(r ) with binding energy Eb, which is equal to 
the one introduced in Sec. [TTJ and the COM wave func- 
tion ^com{R) that experiences the kinetic and potential 
energy of a particle of mass 2m, T-L 2m {R ) ■ 

Consequently, two atoms in an OL at a Feshbach res- 
onance are described by the coupled equations 



H m (f 1 )+H m (f 2 ) + F(r 
W(r) 

Mr)^COM(R) 



W(r) 
U 2m (R) + E b 

^6(»0*COMCB) 



= E 



(20) 



As is usually done for Hubbard models the Hamilto- 
nian is reformulated in the basis of Wannier functions 
of the OL. However, in order to include effects of higher 
Bloch bands and their couplings due to the presence of 
the RBS the basis is not restricted to the first Bloch band. 
In the following the simplification of strong transversal 
confinement is considered, i.e. the ultracold atoms only 
occupy the ground state of transversal motion. Let a\ n 
(tlj,n) be the creation (annihilation) operator of an atom 
with Wannier function = Wi inj o,o and b in (&i, n ) the 
creation (annihilation) operator of the RBS with COM 



Wannier function w; 



i,n,0,0 



The Hamiltonian 



in second quantization that is equivalent to the coupled 
equations (|20|) expanded in the Wannier basis is given as 



(wi, n Wj, m \V \w k ,pWi, q ) a\ n (£ m ak,pai, q 



H = E E (Wi,n\Mm \u>j,m) a\,n a j,r. 



1,3 



\ E E 

i,j,k,l n,m,p,q 
E E (v*i,nl %m \Wj, m ) + E^j b\ n b hm 

^ E E ( W i,n w J,m\ W \w Kp ip b ) (a\ n a] m b k , p + h, 

(21) 



V2 



i,j,k n,m,p 



Note the factor 1 / y/2 before the atom-molecule coupling, 
which has to be included to ensure that the matrix ele- 
ments of the Hamiltonian are equal in first and second 
quantization |31) . 



We want to emphasize that Eq. ( 20 ) and thus the sec- 



ond quantized Hamiltonian ( 21 1 are only valid if not more 



than two atoms interact. For more atoms important ef- 
fects such as losses or the appearance of Efimov states 
cannot be correctly reproduced. 

The following simplifications and approximations are 
introduced: 

1. The Hamiltonians Ji rn and iH. 2m do not cou- 
ple different Bloch bands, since the Bloch func- 
tions Wi tn and Wi tn are eigenstates of T-L m and 
H2 m , respectively. For example, for H m holds 

(Wi,n\ Hm \w 3 ,k) = (w it n\H m \Wj >n ) 5 nk . 

2. Only next-neighbor coupling is considered, i.e. 

E E ( W i,n\H m \w j>n ) al^Gj.r, 



1,3 n 



m | "^,771/ 



E b ) 4,n b 3,™ 



i,j n,m 

i n (i.j) n 

+ E E ( f « + b lnKn "EE fAnhn 

i n (ij) n 

where (• • • ) denotes summation over nearest- 
neighbor lattice sites, e„ = {w\ n \'hL m \wx n ), £ n = 



{Wl,n\ %2ro \wi, n ), Jn 
Jn = - (Wl,n\ U2?n \w 2 . 



(wiJiH m \w2,n), and 



3. The interaction potential is replaced by the Fermi- 
Huang pseudo potential V(r) —> 4lTh ^ b e §(j?jd r 
that reproduces the same background scattering 
length abg as the full open-channel interaction po- 
tential. For small background scattering length 
only onsite-interaction is taken into account, i.e. 

E E ( W i,nWj, m \V\Wk,pWl, q ) al,n a j,m a fc,P a i,q 
i,j,k,l n,m,p,q 



E 



i n,m,p,q 



with f7, 



n,m,p,q 



= (wi, n Wi t m\ V \Wl,pWl,q) = 



47vh flb 



/ dx dy dz w ^ n w ,m wa, P w , q - 



4. The coupling to the molecule happens only at short 
distances, i.e. on the length scale of the lattice and 
the transverse harmonic confinement, thus one can 
replace 



W(r)ip b (r) -t gS(r) 



(22) 



() 



where the coupling strength g has to be adapted to 
match the behavior of the system under consider- 
ation. Including only next-neighbor coupling leads 
to the simplification 



^ ^ (u>i )7 it0,-,z|W ipb) (<4,n a ],i b k,p + h.c^ 

E E 9 { :; n k r k) ( a ln<4,i b k, P + h - c ) > 



i,j,k n,l,p 



{ij,k) n,l,p 



with 



9n'il = 9 I dxdydz w,. n w jtl w 0iP . 



(23) 



Due to the symmetry of the Wannier functions the 
onsite coupling obeys the selection rule 



(0,0) 



= for n + I + p even. 



Employing the above simplifications and approxima- 
tions the BH Hamiltonian reduced to the first N Bloch 
bands is given as 

N N 

H B H = X] X] £ " a i,n a i,n ~ E E J " a i,n a j,n 
i n=l (ij) n=l 

1 N t t 

i n,l,p,q—l 
N N 

+ E E ( f » + E ») b Uhn "EE^ 6 l A" 
i n=l (ij) n=l 

1 N 



(ij,k) n,l,p—l 



(24) 



easily account for the altered resonance parameters by in- 
troducing an effective coupling strength and an effective 
bound-state energy, 



9 -> 3cff 



Uirh 2 a hg AE 



2/t 



S/, -> Bi 



b.cff 



-Ere 



(25) 
(26) 



that lead to the correct resonance width Ai? and reso- 
nance energy E les . In the following the index "eff" will 
be suppressed keeping however in mind that g and Eb 
are not equivalent to the physical coupling strength and 
the physical energy of the RBS. 




IV. PROBLEM OF REPRESENTING A 
DELTA-LIKE COUPLING WITHIN THE 
BOSE-HUBBARD MODEL 

The coupling of the open channel to the bound state 
as described by Eq. ( 22 1 seems to be a crude approxi- 



mation. Indeed, as discussed in Sec. [H| a more general 
form of a short-range coupling to the bound state is of 
the form W(f)tpb(r) = \f^n^ir — a*)S(f). While one 
can associate g with ^/iirja* the coupling \/4^T^rS(r) 
automatically vanishes for the chosen single-atom basis 
states. In fact it vanishes for any basis that conforms 
to a scattering length a — 0. Hence, the presented BH 
model can only conform to a FR with 7 = 0, a* — > 00, 
and 7a* = const. This is only the case for qq = a [see 



Eq. (15)] and results according to Eqs. (13) and (14) in 

2fJ, ha') 2 



a resonance width AE 



-^r — and a res- 



h 2 abg 47r/i 2 abg 

onance detuning SE = 0. For FRs with 7^0 one can 



FIG. 2. (Color online). Energy spectrum as a function of the 
resonance energy E rcs for r\ = 4, ab g = 0.04ah o and from top 
to bottom AE = (1,4, 16)ftoj. This results in the coupling 
energies x [ see Eqs. ( 22 1 and ( 27 1] given in the graphs. The 
analytic eigenenergies (dots) obtained by Eq. ( |16[ ) are com- 
pared to the eigenenergies of the BH model (lines) for two 
Bloch bands (left column) and four Bloch bands (right col- 
umn) included. Including only two Bloch bands (left column) 
the analytic eigenstate with energy « 12huu is not reproduced 
by the BH model. 

In Fig. [2] the energy spectra in an anisotropic harmonic 
trap of several FRs of different widths are compared to 
the corresponding result of the effective BH model. The 
trapping frequencies are uj y = u) z = r]U), with 77 = 4 and 
u! the trapping frequency in x direction. In the harmonic 
trap the Wannier functions of the BH model are replaced 
by harmonic-oscillator eigenfunctions. On the left side 
two Bloch bands are included and the RBS appears in 
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two different COM states while the unbound atoms can 
occupy three different trap states [(i) both atoms in the 
first band at 9hu>, (ii) one atom in the first and one in the 
second band at lOHw and (iii) two atoms in the second 
band On the right side four Bloch bands are in- 

cluded with correspondingly more molecular states and 
trap states. 

As a measure for the coupling strength the energy 



X 



(o,o) 
"1,1,1 



(27) 



is introduced [see Eq. (23)]. The avoided crossing be- 
tween the lowest bound state and the first trap state has 
a splitting energy of rj 2\- 

For a relatively narrow FR with an effective coupling 
strength x — OAHuj the agreement between the BH model 
and the analytic result is very good independently of the 
number of Bloch bands included. For the broader FRs 
with x — O.&Huj and x = one can make two obser- 

vations: (i) Trap states (i.e. states above the bound state 
threshold of 9huj) quickly approach to the analytic results 
for an increasing number of Bloch bands, (ii) The dis- 
agreement between analytic and BH results of the bound 
states does not decrease with the number of Bloch bands. 

Obviously, the variational principle does not hold for 
the bound state as an insufficient basis leads to an energy 
lower than the correct bound state energy. Moreover, 
by increasing the basis the already incorrect bound-state 
energy becomes even lower and the disagreement to the 
correct result increases. Though less severe, the same 
problem also appears for trap states. For example, the 
first trap state in the last row in Fig. [2] lies below the 
correct energy if four Bloch bands are included. 




radial wavefunction (j> E (r) 

expansion to 3 rd order 

expansion to 10 th order 



1 2 3 4 5 

r [units of a h0 ] 

FIG. 3. (Color online). The radial wave function 4>E{r) 
in a spherical harmonic trap introduced in Eq. |5| is com- 
pared for the REL energy E = 2.5hui to its expansion 
<?W = Yj1=o( ( t >n \ c t 3 E)4>n{r) to different orders N, where <j> n 
is the radial wave function of the non-interacting system with 
radial momentum I = and energy (2n + %)tkj. Since all 
non-interacting radial basis functions are zero for r = the 
expansion cannot reproduce the behavior of 4>b{t) for r — > 0. 
This is important, since the coupling to the bound state is 
proportional to (/>'e(0) or 0cx P (O), respectively. 

The reason for this insufficiency of the basis to conform 
to the behavior of a delta-like coupling is related to the 



problem of a missing coupling of the form ^/A^"fr8{r): 
the two-particle basis states are a = wave functions. 
However, a = basis functions can represent the full 
wave function only for r > but not for r — > (see 
Fig. p$J. While for ordinary interaction potentials the 
value of the wave function at r = is irrelevant, for 
zero-range potentials it is decisive. The problem is es- 
pecially severe for the open-channel bound state, which 
appears for positive scattering lengths. For E — » —00 
one has |</>.e(0)| oc (— E) 1 ^ making its representation by 
a = basis functions for decreasing energy more and 
more problematic. 

For weak coupling the problem is less severe as eigen- 
states that differ significantly from the background trap 
states are predominantly bound states with different 
COM excitations, which are well reproduced by the BH 
model. For strong coupling, however, the bound state is 
admixed to many states in the spectrum (see Sec. |n| . 
Since the bound-state admixture for a certain eigen- 
state is thus lower, a good representation of the open- 
channel wavefunction is important also for large scatter- 
ing lengths. 

The described problem does not only arise when using 
non-interacting a = basis states. For any finite expan- 
sion of the radial wave function ^ exp (r) = ^2 c n <j) n (r) in a 
superposition of basis functions with a specific scattering 
length at, [i.e. a& = — 4> n (0) / 4>' n (0)} the scattering length 
of the expansion yields 



^xp(°) 



Ob . 



(28) 



Hence, the wave function (f> exp (r) cannot adapt to a 
change of the scattering length induced by a short-range 
coupling. Especially, since the scattering length at a FR 
is energy dependent these expansions cannot reproduce 
the correct eigenenergies and eigenstates. 



V. DRESSING OF COUPLING STRENGTH 
AND BOUND-STATE ENERGY 

To circumvent the problem of the wrong representa- 
tion of a zero-range coupling one can replace it by a 
finite-range coupling. To this end one usually considers 
the Fourier transform of the problem and regularizes the 
delta-like interaction by introducing a high-momentum 
cut-off A. Thereupon the coupling parameter is renor- 
malized [32]. Taking the limit A — > 00 the finite-range 
coupling converges towards a zero-range coupling. How- 
ever, for an interaction with a range of d/N where d is 
the lattice spacing more than N Bloch bands have to be 
included to converge the energies [7]. 

Here we want to take a different approach with no need 
to include more Bloch bands to reproduce the correct 
bound-state energies. Provided with the analytic solu- 
tion in the harmonic trap a dressed bound state is intro- 
duced, which reproduces the correct energy spectrum in 
the harmonic trap at least in the important energy range 



of the first Bloch band. We use the fact that the full 
bound state (the combination of the closed-channel and 
open-channel bound state) falls off rapidly for increas- 
ing internuclear separation. Hence, the bound state does 
hardly probe the anharmonic parts of the potential and 
the dressed bound state can be equally used for (anhar- 
monic) OLs. 




£ res [units of hoj] 

FIG. 4. (Color online). Results of the dressed BH model for 
four Bloch bands with tibg = 0.04ai lo and AE = 16ftu. Top: 
Dressed bound-state energies and dressed coupling strength 
as a function of -E rcs . Bottom: Comparison of the ana- 
lytic energy spectrum (dots) with the energies of the dressed 
BH model (solid lines) and the undressed BH model (dotted 
lines). 

More concretely, the dressed bound state is introduced 
in the following way: The RBS in the first band (for 
which the COM wave function is a Wannier function 
of the first band) couples predominantly to two atoms 
in the first band leading to the lowest avoided cross- 
ing in the spectrum. The two corresponding eigenen- 
ergies are given by a sum of the lowest COM energy 
E COM [#com = ^ n _ i + 77)] and the two lowest so- 
lutions Ei , E2 of the REL motion eigenenergy relation 
( 16 1 which depend on the bound-state energy Ef, = E ICS . 



In order to match the energies of this avoided crossing 
the bound-state energy E^ and the coupling strength g 
are replaced by dressed parameters E^ — ► Ej(E Ies ) and 
9 9d- The two parameters are determined by a least 
square fit to the energies Ei + E$ OM and E 2 + Efi OM . 

To match the energies E x + E% OM with n = 2,3,... 
of bound states in higher Bloch bands, dressed bound- 
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state energies E d 2) (E Ies ), E { d v (E rca ) 
which are also determined by a least square fit. The up- 



p(3), 



are introduced, 



FIG. 5. (Color online). Energy spectrum of two atoms in an 
OL with lattice depth Vl = 5E r — l.lhtu consisting of three 
lattice sites with periodic boundary conditions. Excitations 
in transversal direction are frozen out by choosing transver- 
sal trapping frequencies u) v = uj z = 3.8cj, where ui is the 
frequency of the harmonic approximation of a lattice site in 
x direction. The resonance parameters are ab g = 85 a.u. = 
9.0T0 _3 d, and AE = 2A.2frw which corresponds to a coupling 
strength of \ = 1-66 fiw = 1.48 Vl (See also the right graph 
in Fig. [9] with the same lattice parameters and resonance pa- 
rameters). The comparison of the eigenenergies of the dressed 
BH model (solid lines) and the undressed BH model (dashed 
lines) each with four Bloch bands included, shows that again 
both models disagree especially for the bound states, while 
the differences for the trap states are small. The inset shows 
a magnification of the spectrum close to the crossing of ex- 
cited bound states with the lowest Bloch band. 



per branches of the avoided crossings with bound states 
in higher Bloch bands lay above the first Bloch band. 
Therefore, their correct representation is less relevant and 
we do not need to introduce also band-dependent dressed 
coupling strengths. 

In Fig. Q the dressed energies E d l) , E [ d ] , Ef ] , E { d i} and 
gd and the corresponding corrected spectrum are shown 
for the four-band BH model with a^ g = 0.04ai lo and 
AE = IQhio (same parameters as for right bottom graph 
in Fig. [2j. Evidently, the dressing of the bound states 
becomes relevant for a resonance energy E Tes < 5Hto, 
but is already visible for E ICS < 10 hiu . Since only a 
band-independent dressed coupling strength was intro- 
duced, the repulsive branches above the first Bloch band 
with an energy above 10Hu> are not fitted to the exact 
results. Correspondingly, slight deviations between the 
exact energies and the dressed BH energies appear for 
these states, while the first repulsive branch is correctly 
reproduced. 

The introduced dressed parameters can now be used 
to determine the energy spectrum of two atoms in an 
OL. In Fig. [5] the spectrum of the dressed and undressed 
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BH model of two atoms in a small OL consisting of 
three lattice sites are compared for a coupling energy 
of x — 1.66 Tied = 1.48 Vl- In contrast to the purely 
harmonic trap, the energies of the bound states and the 
trap states split due to tunneling. If the molecular bound 
states are not in resonance, i.e. for E TCS < 0, the trap- 
state energies form bands of increasing widths around 
8Ahw, 9.1huj, 9.8Hlj, and IOAHuj. For resonance ener- 
gies -E res > the bound states cross with the trap states 
leading to a plethora of avoided crossings. In the ul- 
tracold regime especially the avoided crossings with the 
first band are of relevance. These appear due to the 
next-neighbor coupling of the molecular state with the 
atomic states [33]. As shown in the inset of Fig. [5] the 
width of these avoided crossings decreases with the COM 
excitation energy of the RBS. The comparison between 
the dressed and the undressed BH model shows that also 
in the OL the energies disagree especially for the bound 
states, while the energy differences for the trap states are 
small. 



VI. NON-PERTURBATIVE DETERMINATION 
OF STATIONARY AND DYNAMICAL STATES 

In the following the results of the BH model shall be 
compared to non-perturbative calculations for two atoms 
at a FR in an OL consisting of two lattice sites. In order 
to do so an approach described in |21j is used, which al- 
lows for finding the stationary solutions of the two-body 
problem with arbitrary isotropic single- channel interac- 
tion potentials. On the basis of the stationary solutions 
the method described in [S3] is used to determine the 
timc-dependent wavefunction during a perturbation of 
the lattice potential. 

Since the lattice potential couples REL and COM mo- 
tion and the interaction couples the motion in x, y and z 
direction all six coordinates of the problem are coupled. 
An extension to the coupling to an additional channel 
describing the COM and REL motion of the molecular 
bound state would make the solution very cumbersome. 
Instead, the freedom of the choice of the interaction po- 
tential is used to realistically mimic a two-channel prob- 
lem by a square-well interaction potential. The potential 
supports bound states that are coupled by a barrier to 
the scattering states. In the following it is shown that 
this potential leads to an energy dependence of the scat- 
tering length, which is in very good agreement to the one 
of a two-channel description in Eq. (12). This is already 



sufficient to realistically mimic a FR since, as shown in 
Sec|H] the energy dependence of the scattering length 
fully determines the energy spectrum. 
The square-well potential is defined as 



-Vq for r < t*o 
V(r) = < +Vi for r < r < n 
elsewhere 



(29) 




FIG. 6. (Color online). Sketch of the square well potential 
(thick blue) . The resonant character of the potential is due to 
the coupling of a bound state with energy Eb (gray dashed) 
to an unbound state with energy E (gray dashed). The ac- 
cording wave functions are sketched by black thin lines. For 
Eb ~ E the scattering length changes resonantly. 



been used to study effects of the energy-dependence of 
the scattering length on the BEC-BCS crossover [33]. For 
sufficiently large Vo the potential supports a bound state 
behind a potential barrier of height V\ and width r± — r$. 
An atom pair that collides with an energy E = h 2 k 2 / (2/i) 
scatters resonantly, if E is close to the bound-state en- 
ergy. 

Introducing dimensionless variables p = r/ri, d = 
ra/ri, k = kr\, v = Vq/E q , and v\ — Vi/E with 

Eq = tP^i the solution of the Schrodinger equation for 

E > is given as 



C sin(fco/o) 
4>{p) = \ Ae klP + Be' klP 
sin [up + ip(k)) 



for p < d 
for d < p < 1 
elsewhere, 



(30) 



with fco = \/vq + k 2 and k\ — \/vi — k 2 . 

In the case of pure s-wave scattering one has K <C 1 
so that one can make, e.g., the replacements sin(/c) — > k 
and cos(k) — > 1. Eliminating A, B, and C by demanding 
that the wavefunction is continuous and differentiable the 
scattering length can be obtained as 



o(k 2 ) _ taxi(p 1 + ' 



Ere 



with 



ce — ft 

= e 2dkl [ho cos(dk ) - ki sin(dfco)] ■ 
= e 2kl [ko cos(dk ) + k\ sm(dk )} . 



(31) 

(32) 

(33) 
(34) 



From the functional behavior of Eq. ( |31[ ) one can deter- 
mine the corresponding parameters of the FR, i.e. E xes , 
AE, and ab g - The resonance positions of a{n 2 ) are given 



by the roots of K 



be called k 2 b 
evaluates to 



e ( K rcs) 



3 (k 2 ). The smallest root shall 
Hence, the resonance position 



with Vb,Vi > (see Fig. |6|. This potential has also 



E r . 



E k 2 



(35) 
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According to Eq. ( 12 ) the scattering length is zero if E 
E, 



AE. Be kq the solution of 1 + e res (Ko) 
closest to Kj-es then 

AE = E (k 2 cs -k 2 ). 



length one could add another square well with V < in 
that is front of the potential in Eq. ( 29 ) . However, here the focus 



(36) 



In order to determine the value of the background scat- 
tering length dbg, e rC s is expanded linearly in k 2 around 
the resonance position, yielding 



,(« 2 ) 



with 5 



k 2 



d{K 2 



S(k 2 



(37) 
(38) 




FIG. 7. (Color online) Energy-dependent scattering length of 
the square-well potential (dots) and approximation according 
to Feshbach theory (thin red) for Vi = 70 and ro/n = 0.6. 

For k — > K r0 s the scattering length evaluates according 



to Eq. (31) and Eq. (37) to 



a(n 2 ) 

n 



1 ( K 2 
5-1 \^res 



By comparing with the behavior of Eq. ( 12 ) for E 

a = ai, g AE/(E rcs — E), one finds 



E„ 



a hf . 



Eq 



AE{5-1) 



(39) 



E rQ 



(40) 



For non-resonant background scattering the wavefunc- 
tion simply falls off exponentially for r < r±. Therefore 
fl bg J; ri. Since the potential mimics an s-wave reso- 
nance, the choice for r± is limited to fcrj <C 1 and for 
energies E ~ fko to r% <C a\ lQl allowing only for rather 
small positive background scattering lengths. On the 
other hand, one can freely choose -E res and AE by an 
appropriate choice of the parameters Vq and V\ , respec- 
tively. In order to also control the background scattering 



lies on the coupling to the RBS and not on the value of 

«bg- 

In Fig. |7ja(K 2 ) is shown for an exemplary square-well 
potential with d — 0.6 and V\ — 70. The values of a(n 2 ) 



according to Eq. (31) and its approximation 



o bg 1 - 



AE 



E te 



E 



(41) 



with the parameters according to the equations ( 35 1 , 



(36), and (40) agree almost perfectly, showing that the 



square-well potential reproduces very well the behavior 
of a FR. 



VII. COMPARISON OF BOSE-HUBBARD 
MODEL TO NON-PERTURBATIVE 
CALCULATIONS 

A. Energy spectrum 

Equipped with the possibility to model FRs with a 
single- channel potential we can apply the ab-initio ap- 
proach introduced in |21j to determine the energy spec- 
trum of two atoms at an FR in a small OL with a lattice 
spacing of d = 500 nm. Within the numerical approach 
one can expand the OL potential in all directions to some 
arbitrary order. Again, to avoid unnecessary complexity 
the OL is expanded to harmonic order around y = z = 
in y and z direction and to 12-th order around x = n/2 
in x direction. The lattice depth in y and z direction 



is chosen sufficiently large (u> y 



0J Z 



3.8w where to is 



the trap frequency of the harmonic approximation of the 
lattice wells in x direction) such that excitations in these 
directions can be ignored. The resulting double-well po- 
tential in x direction is shown in Fig. [8j For large 
lattice depths the spectrum converges to the one of two 
uncoupled harmonic traps. In order to probe the accu- 
racy of the BH model a relatively small lattice depth of 
Vl = 5E r = l.lhu) is chosen in x direction. For this low 
lattice depth excited states in higher Bloch bands probe 
parts of the potential that significantly deviate from an 
ordinary lattice potential Vl sin 2 (k x). Therefore, the 
correct single-atom states deviate significantly from or- 
dinary Wannier functions. This insufficiency can be cor- 
rected for by replacing the ordinary Wannier basis by 
a basis constructed from single-atom eigenstates in the 
double well. For each band p the left and right Wan- 
nier functions are constructed by superpositions of the 
ri-th symmetric eigenstate with energy £^p 6ven ^ and the 

n-th anti-symmetric eigenstate with energy £ , 4° dd ' ) - The 
corresponding atomic Wannier functions of the first four 
Bloch bands are shown in Fig. [Hj As one can see they 
are neither symmetric nor anti-symmetric so that any se- 
lection rule for the BH parameters (such as that of the 
coupling between the open and the closed channel) of 
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FIG. 8. (Color online) Double-well potential (thick, solid) 
used in the ab initio calculations and corresponding full lat- 
tice potential VL sin 2 (fco:r) (thick, dashed). The Wannier func- 
tions of the atoms in the BH model are depicted for bands 
one to four (red solid, green dashed, blue dotted and orange 
dot-dashed) alternately for the right and the left well. Al- 
ready above the first band they clearly probe regions where 
the double-well potential significantly differs from the full lat- 
tice potential. Horizontal lines mark the onsite energies of 
bands one to four. 



harmonic counterparts (see Fig. [8]). Nevertheless, the 
dressed bound-state energies and the dressed coupling 
strength lead to a drastic improvement of the undressed 
results in all three cases shown in Fig. [9j In general, the 
dressed parameters should lead to an improvement as 
long as the couplings of the bound states to trap states 
that probe anharmonic parts of the potential, i.e. with 
energies above E = Vl, is negligible. Approximately, 
for x > Vl this is not the case any more since at the 
avoided crossing of the lowest bound state with the low- 
est trap state an energy regime above Vl is entered. In- 
deed, considering the spectrum with the largest coupling 
energy \ — 1.48 Vj, = 1.66/tw, the lowest bound-state 
energy of the BH model is slightly lower than that of 
the ab-initio calculations. But still the disagreement is 
surprisingly small. As one can expect the correction of 
the bound-state energies in the third and fourth Bloch 
band is less accurate than that of the first and second 
Bloch band. Already for the lower coupling energies of 
X = 0.36V L = 0.41/ia; and X = 0.67VL = 0.75/iw small 
disagreements between the corresponding eigenenergies 
of the ab-initio calculations and the corrected BH model 
appear. 



the OL does not apply. The onsite energies are given as 

and the hopping parameters as 
. Furthermore, to be sure that 



n(odd) 



^(even)^ 
- En V ) 



t _ l/p(° dd ) 

all errors are solely due to deficiencies of the represen- 
tation of the Feshbach resonance in the BH model also 
next-neighbor (background) interaction is included. 

In Fig. [9] the spectrum of the ab-initio calculation for 
three different coupling strengths is compared to the cor- 
responding dressed and non-dressed BH spectrum. In 
contrast to Fig. [5] the trap states do not appear in en- 
ergy bands due to the reduced size of the system. The 
bound states appear as duplets with one symmetric and 
one antisymmetric COM excitation in x direction. Again, 
excited bound states in higher Bloch bands are able to 
couple to the first trap state (lowest horizontal line) by 
next-neighbor coupling, i.e. the bound state couples to a 
state of one atom in the same well and one in the neigh- 
boring well. For symmetry reasons only the lower bound 
state of each dublet can couple to the lowest symmetric 
trap state [33] , 



Fig. 10 shows a detailed view onto two of these avoided 
crossings around E — 8Ahui for a resonance energy of 
E xes — 2.9hoj and E les = 3.9fku. Clearly, the splitting 
of the avoided crossing and hence also the next-neighbor 
coupling strength is well reproduced by the dressed BH 
model. 

Given the large degree of anharmonicity of the lat- 
tice potential the agreement between the ab-initio spec- 
tra and BH spectra in Figs. [9] and [10] is very good. 
The dressed bound-state energies are obtained from a 
harmonic approximation of the two lattice sites. Al- 
ready in the second Bloch band the potential and there- 
fore the states and energies differ significantly from their 



The coupling of the two atoms in the lowest Bloch 
band to the bound state in the lowest Bloch band leads 
to the appearance of both attractively and repulsively 
interacting states. The energy of the repulsively inter- 
acting state is marked by the red shading in Figs. [9] and 
|10| As one can see for larger and larger coupling energy x 
this state is strongly influenced by bound states in higher 
and higher Bloch bands. If this energy range shall be cor- 
rectly reproduced this sets a lower limit for the number of 
Bloch bands that must be included in the BH model. In 



Fig. 10 one can see that the dressed BH model reproduces 
correctly the energy of the repulsively interacting state 
while the undressed model underestimates its energy. 

As discussed above, the dressed BH model reproduces 
accurately the correct eigenenergies up to coupling en- 
ergies x ~ Vl- This corresponds usually to small up to 
medium FRs. As discussed in Sec. |H] a FR in a har- 
monic trap is broad if abgAi? 3> atohiij. Since x is a 
measure for the energy splitting of the avoided crossing 
of the lowest bound state with the first band, it is com- 
parable to yj ahgAE j (a\ 10 hu>)huj in the harmonic trap. 
Therefore, a broad resonance requires [x/(/iw)] 2 ^> 1. 
Since the BH model is valid for X ~ Vl it can only ac- 
curately describe broad FRs in a very deep lattice with 
(V L /(M) 2 = V L /(4E r ) » 1. 

However, for broad FRs all details of the interaction 
apart from the value of the scattering length for E — > 
are irrelevant (see Sec. |TT[ ) . In this situation there is 
not required to explicitly include the bound state in the 
BH model. Instead, corrected BH models like the one 
introduced in [5] already provide accurate results. 
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FIG. 9. (Color online) Spectra of the ab-initio calculations (dots) and the BH model with usage of the dressed bound-state 
energies and coupling (lines). Also shown are the energies of the undressed BH model (dotted lines). The ab-initio calculations 
include the representation of bound states with many COM excitations. Not all of these bound states are present in the BH 
model that only includes four Bloch bands. For example, in the right graph all ab-initio energies for 7Ahw < E < 8Ahu and 
Eres < are not covered by the BH model. From left to right the parameters ab g = (88, 87, 85) a.u. = (9.3, 9.2, 9.0) • 10 _3 d, 
AE = (1.4, 4.9, 24.2)fiw are chosen. This corresponds to a coupling strength of x = (0-41, 0.75, 1.66) huo = (0.36, 0.67, 1.48) Vl. 
The red shading marks the energy of the repulsively interacting atoms within a single-band approximation. From left to right 
the energy of this state is significantly influenced by the bound state in the second, third and fourth Bloch band demonstrating 
that for stronger coupling bound states in more Bloch bands have to be included to obtain accurate eigenenergies. 
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FIG. 10. (Color online) Zoom on the resonance of the bound 
state in the second Bloch band (right circle) and third Bloch 
band (left circle) with the state of two separated atoms in the 
ground state for \ — 0.75&J. The splitting energies of the left 
resonance (0.04foj) and that of the right resonance (0.06fiu) 
are well reproduced by the dressed BH models. 



B. Time-dependent manipulations 

In the following it is studied how well the BH model 
can predict the dynamic behavior of the system under 
the influence of some time- dependent perturbation 

M 

Vpert(t) = ^2 Vpeitin )/(*), 
i=l 

which acts on each of M identical atoms in the same way. 



Normally, any external potential v pei t(fi) is approxi- 
mately constant on the length scale of the bound state. 
Hence, the perturbation cannot couple the orthogonal 
closed and open channel states. The matrix elements of 
the perturbation of the closed channel evaluate to 

{ipbWj,n\ ^ Vpectjfi ) \wj, m ipb) = dR dr \ip b (r)\ 2 x 



x w it n(R) 



Vpeit I R + 



dRu>i, n (R) J dr \Mr)\ 2 2v pelt (R) (R) 

= 2 (W it n\ Vpert \Wj,m} 

Hence, in second quantization the perturbation is ex- 
pressed as 



Vpert (t) = f(t) f ( Wi - n \ ^P crt \ W i> m "> a ln a 3- 

^ (i j") n,m=l 

N 

+ 2 Y (wi,n\% e n\wj, m )bl n b jtm 

(id) n,m=l 

As usual, only next-neighbor coupling and on-site cou- 
pling are considered and the basis is restricted to the 
first TV Bloch bands. 

In the following, the case of a linear perturbation with 
increasing strength, 

JV 



Vpert (t) — y jXjXt , 
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FIG. 11. (Color online) Sketch of the dynamical behavior 
while accelerating (inclining) the double-well, a) The initial 
state consists of separated atoms (red disks) in the ground 
state of the left and right well. The four molecular states 
in the COM ground state (blue double disk below red disks) 
and in the first excited COM state (blue double disk above red 
disks) are not in resonance, b) Upon inclining the potential 
the energy of an excited molecular state in the left well (dark 
blue) comes in resonance with the energy of the separated 
atoms. The molecular state is occupied and the COM of 
the system moves to the left, c) After a further inclination 
the energy of the excited molecule in the left well comes into 
resonance with the ground-state molecule in the right well. 
By occupying this state the COM of the system moves to the 
right, d) Finally, the molecule on the right well comes into 
resonance with the initial state of two separated atoms and 
the COM of the system moves again to the left. 



is considered, which corresponds to an increasing accel- 
eration of the lattice [35]. The dynamical behavior due 
to Vpert is governed mainly by two effects: (i) The linear 
perturbation leads to a coupling between Wannier func- 
tions of odd and even symmetry, i.e. between bands with 
odd and even quantum numbers, (ii) The energy of the 
states at each lattice site are shifted proportionally to the 
product Xdj and thus depend on the site number j. 

Of course, the dynamical behavior also strongly de- 
pends on the value of the resonance energy E ros . For 
the dynamical studies a resonance energy is chosen such 
that an inclination leads to the resonant next-neighbor 
coupling of two separated atoms in the ground state to a 
bound state in the first and second Bloch band. The cor- 
responding dynamical behavior is sketched in Fig. |11| As 
one can see the COM movement of the system upon ac- 
celerating the lattice depends crucially on the resonance 
energy, i.e. the energy of the RBS. Depending on the 
bound state and its COM excitation that comes into res- 



onance the system can move against the direction or in 
direction of the acceleration. A precise representation 
of the system is thus necessary to predict the mobility 
behavior of two atoms at a Feshbach resonance. 

Fig. [l2] shows the projections | (n\ ^(t)) | 2 of the time- 
dependent wavefunctions ^(i)) onto the eigenstates |ra) 
of the unperturbed system for a slow inclination with 
A = 0.0003^^. If the perturbation would be suddenly 
switched off, the projections give the probability of find- 
ing the system in the corresponding eigenstate. For the 
same three coupling energies as shown in Fig. [9] the qual- 
itative agreement between the result of the ab-initio ap- 
proach (upper row) and the dressed BH model (middle 
row) is very good. As is visible in Fig. 11 initially the 
bound state in the second Bloch band is slowly occupied. 
After t « 1300h/E r this bound state gets into resonance 
with the bound state in the first Bloch band which is 
then occupied. After t « 1500h/E r the main occupa- 
tion moves back to the initial state. Additionally to the 
behavior described in Fig. [TT] the inclination leads to a 
strong coupling of the bound states in the first and sec- 
ond Bloch bands on each lattice site. Due to the large 
energy separation of these states this coupling leads to 
fast oscillations of the population of the eigenstates. 

In order to examine the quantitative agreement be- 
tween the ab-initio and dressed BH results the time- 
dependent COM motion of the system X 
has been determined. As one can see in the lower row 
in Fig. [TT] the quantitative agreement is very good for 
the smallest coupling energy \ — OAlHu, For the larger 
coupling energies especially the fast oscillations appear- 
ing after t ps 1200h/E r are less accurately reproduced by 
the dressed BH model. The phase shift and altered fre- 
quency of the oscillations is mainly due to a small under- 
estimation by about 1% of the coupling strength between 
the stationary eigenstates within the dressed BH model. 
In contrast to the dressed BH model, the undressed BH 
model leads even for small coupling energies to a dynam- 
ical behavior significantly disagreeing from the one of the 
ab-initio calculations. 



VIII. CONCLUSION 

We have introduced a Bose-Hubbard model with 
dressed bound states and a dressed coupling strength, 
which can be used to accurately determine the stationary 
and dynamical wavefunctions of two atoms in an optical 
lattice at a Feshbach resonance. The dressed parame- 
ters, which can be straightforwardly obtained from the 
analytically known solution of a Feshbach resonance in a 
harmonic trap, allow one to obtain an accurate solution 
with including only a small number of Bloch bands. The 
dressing avoids the problem that the eigenenergies, ob- 
tained by a finite expansion of the solution in single-atom 
basis states, do not converge to the correct eigenenergies 
in the presence of a delta-like coupling to the bound state. 
Hence, the introduced method permits to determine ac- 
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FIG. 12. (Color online) Dynamic behavior of two separated atoms in the ground state of the double-well potential during an 
inclination of the lattice for different coupling energies x an d resonance energies E Tas . At t = i cll d = 2000h~/E r each atom 
experiences a perturbation of W = Q.lhujx/d which suffice to bring both the first and second bound state into resonance (see 
Fig. 1 1 1 [ ) . The projection of the time-dependent wave function onto the eigenstates |n) of the unperturbed system is 

shown in the first row [ab-initio results) and the second row (results of the dressed BH model) using the same color coding as 
in Figs. [9] and 
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In the lowest row the mean COM position (^(t)! X \^f(t)) is shown for the ab initio results and the dressed 
and undressed BH model. The insets show a magnified region of the beginning of the fast oscillations between 1200h/E r and 
U00h/E r . 



curate solutions without a regularization of the potential 
and a numerically demanding expansion of the solution, 
e.g., in Bloch functions or Wannier functions of many 
Bloch bands. The solution of this problem should be rel- 
evant to many approaches that seek to describe strongly 
interacting atoms via a multi-band Hubbard model. 

Comparisons to eigenenergies and time-dependent 
wavefunctions obtained from a non-perturbative ap- 
proach have shown that the method is accurate as long 
as the coupling energy is smaller or comparable to the 
lattice depth. Furthermore, we have described a possi- 
bility to realistically mimic FRs within non-perturbative 
single-channel approaches by using a square- well interac- 
tion potential. 



We believe that the approach is applicable not only 
to optical lattices but to various kinds of anharmonic 
trapping potentials. The introduced methods should be 
therefore a valuable tools for investigating the exciting 
physics of Feshbach-interacting atoms in various poten- 
tials and to interpret corresponding experimental find- 
ings. 
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